7-6 shortness of breath

Shortness of breath.

Authors

Affiliations

James Totterdell

University of Sydney

Rob Mahar

University of Melbourne

Load packages

library(ASCOTr)
library(tidyverse)
library(lubridate)
library(kableExtra)
library(patchwork)
library(cmdstanr)
library(posterior)
library(bayesplot)
library(ggdist)
library(lme4)
library(broom)
library(broom.mixed)
library(bayestestR)

theme_set(theme_classic(base_size = 10, base_family = "Palatino") +
  theme(panel.grid = element_blank(),
        strip.background = element_blank()))

bayesplot_theme_set(theme_set(theme_classic(base_size = 10, base_family = "Palatino") +
  theme(panel.grid = element_blank(),
        strip.background = element_blank())))

color_scheme_set("red")
options(digits = 4)

Prepare datasets

all_dat <- readRDS(file.path(ASCOT_DATA, "all_data.rds"))

# FAS-ITT
fas_itt_dat <- ASCOTr:::make_fas_itt_set(all_dat)
fas_itt_nona_dat <- fas_itt_dat |>
  filter(!is.na(out_sob))

# ACS-ITT
acs_itt_dat <- ASCOTr:::make_acs_itt_set(all_dat)
acs_itt_nona_dat <- acs_itt_dat |>
  filter(!is.na(out_sob))

# AVS-ITT
avs_itt_dat <- ASCOTr:::make_avs_itt_set(all_dat)
avs_itt_nona_dat <- avs_itt_dat |>
  filter(!is.na(out_sob))

Load models

logistic_mod <- compile_cmdstanr_mod(
  file.path("binary", "logistic"), dir = "stan")
logistic_site_epoch <- compile_cmdstanr_mod(
  file.path("binary", "logistic_site_epoch"), dir = "stan")

Code

make_summary_table_anticoagulation <- function(dat, format = "html") {
  tab <- dat %>%
    mutate(CAssignment = factor(
      CAssignment, 
      levels = c("C1", "C2", "C3", "C4"),
      labels = str_replace(intervention_labels()$CAssignment[-1], "<br>", " "))) %>%
    group_by(CAssignment) %>%
    summarise(
      Patients = n(),
      Known = sprintf(
        "%i (%.1f)", sum(!is.na(out_sob)), 100 * mean(!is.na(out_sob))),
      Missing = sprintf(
        "%i (%.1f)", sum(is.na(out_sob)), 100 * mean(is.na(out_sob))),
      `Shortness of breath day 28` = sprintf(
        "%i (%.1f)", sum(out_sob, na.rm = TRUE), 100 * mean(out_sob, na.rm = TRUE)),
    ) %>%
    bind_rows(
      dat  %>%
        group_by(CAssignment = "Overall") %>%
        summarise(
          Patients = n(),
          Known = sprintf(
            "%i (%.1f)", sum(!is.na(out_sob)), 100 * mean(!is.na(out_sob))),
          Missing = sprintf(
            "%i (%.1f)", sum(is.na(out_sob)), 100 * mean(is.na(out_sob))),
          `Shortness of breath day 28` = sprintf(
            "%i (%.1f)", sum(out_sob, na.rm = TRUE), 100 * mean(out_sob, na.rm = TRUE)),
      )
    ) %>%
    mutate(CAssignment = fct_inorder(CAssignment)) %>%
    gather(key, value, -CAssignment, factor_key = T) %>%
    spread(key, value)
  colnames(tab)[1] <- "n (\\%)"
  if(format == "latex") {
    colnames(tab) <- linebreak(colnames(tab), align = "c", linebreaker = "<br>")
  }
    kable(tab,
      format = format,
      align = "lrrrrr",
      escape = FALSE,
      booktabs = TRUE
    ) %>%
    kable_styling(font_size = 10, latex_options = "HOLD_position")  %>%
    row_spec(nrow(tab), bold = TRUE)
}


make_summary_table_antiviral <- function(dat, format = "html") {
  tab <- dat %>%
    mutate(AAssignment = factor(
      AAssignment, 
      levels = c("A1", "A2"),
      labels = str_replace(intervention_labels()$AAssignment[-1], "<br>", " "))) %>%
    group_by(AAssignment) %>%
    summarise(
      Patients = n(),
      Known = sprintf(
        "%i (%.1f)", sum(!is.na(out_sob)), 100 * mean(!is.na(out_sob))),
      Missing = sprintf(
        "%i (%.1f)", sum(is.na(out_sob)), 100 * mean(is.na(out_sob))),
      `Shortness of breath day 28` = sprintf(
        "%i (%.1f)", sum(out_sob, na.rm = TRUE), 100 * mean(out_sob, na.rm = TRUE)),
    ) %>%
    bind_rows(
      dat  %>%
        group_by(AAssignment = "Overall") %>%
        summarise(
          Patients = n(),
          Known = sprintf(
            "%i (%.1f)", sum(!is.na(out_sob)), 100 * mean(!is.na(out_sob))),
          Missing = sprintf(
            "%i (%.1f)", sum(is.na(out_sob)), 100 * mean(is.na(out_sob))),
          `Shortness of breath day 28` = sprintf(
            "%i (%.1f)", sum(out_sob, na.rm = TRUE), 100 * mean(out_sob, na.rm = TRUE)),
      )
    ) %>%
    mutate(AAssignment = fct_inorder(AAssignment)) %>%
    gather(key, value, -AAssignment, factor_key = T) %>%
    spread(key, value)
  colnames(tab)[1] <- "n (\\%)"
  if(format == "latex") {
    colnames(tab) <- linebreak(colnames(tab), align = "c", linebreaker = "<br>")
  }
    kable(tab,
      format = format,
      align = "lrrrrr",
      escape = FALSE,
      booktabs = TRUE
    ) %>%
    kable_styling(font_size = 10, latex_options = "HOLD_position")  %>%
    row_spec(nrow(tab), bold = TRUE)
}



make_out_sob_table_C <- function(dat, format = "html") {
  tab <- dat %>%
    mutate(CAssignment = factor(
      CAssignment, 
      levels = paste0("C", 0:4),
      labels = intervention_labels()$CAssignment)
    ) %>%
    group_by(CAssignment) %>%
    summarise(
      Randomised = n(),
      `Outcome missing` = sprintf("%i (%.1f)", sum(is.na(out_sob)), 100 * mean(is.na(out_sob))),
      `Outcome observed` = sprintf("%i (%.1f)", sum(!is.na(out_sob)), 100 * mean(!is.na(out_sob))),
      `Died within 28 days` = sprintf("%i (%.1f)", sum(out_sob, na.rm = TRUE), 100 * mean(out_sob, na.rm = TRUE)),
    ) %>%
    bind_rows(
      dat  %>%
        group_by(CAssignment = "Overall") %>%
        summarise(
          Randomised = n(),
          `Outcome missing` = sprintf("%i (%.1f)", sum(is.na(out_sob)), 100 * mean(is.na(out_sob))),
          `Outcome observed` = sprintf("%i (%.1f)", sum(!is.na(out_sob)), 100 * mean(!is.na(out_sob))),
          `Died within 28 days` = sprintf("%i (%.1f)", sum(out_sob, na.rm = TRUE), 100 * mean(out_sob, na.rm = TRUE)),
      )
    ) %>%
    mutate(CAssignment = fct_inorder(CAssignment)) %>%
    gather(key, value, -CAssignment, factor_key = T) %>%
    spread(CAssignment, value)
  colnames(tab)[1] <- "n (\\%)"
  if(format == "latex") {
    colnames(tab) <- linebreak(colnames(tab), align = "c", linebreaker = "<br>")
  }
    kable(tab,
      format = format,
      align = "lrrrrrr",
      escape = FALSE,
      booktabs = TRUE
    ) %>%
    kable_styling(
      font_size = 9, 
      latex_options = "HOLD_position")  
}

Descriptive

Relationship anti-coagulation to outcome

save_tex_table(make_summary_table_anticoagulation(
      acs_itt_dat,
      "latex"), 
      "outcomes/secondary/7-6-anticoagulation-summary")
make_summary_table_anticoagulation(
    acs_itt_dat
)

n (\%)	Patients	Known	Missing	Shortness of breath day 28
Low dose	610	577 (94.6)	33 (5.4)	115 (19.9)
Intermediate dose	613	584 (95.3)	29 (4.7)	110 (18.8)
Low dose with aspirin	283	271 (95.8)	12 (4.2)	59 (21.8)
Therapeutic dose	50	44 (88.0)	6 (12.0)	11 (25.0)
Overall	1556	1476 (94.9)	80 (5.1)	295 (20.0)

Table 1: Primary outcome by anti-coagulation intervention.

Relationship anti-viral to outcome

save_tex_table(make_summary_table_antiviral(
     avs_itt_dat,
      "latex"), 
      "outcomes/secondary/7-6-antiviral-summary")
make_summary_table_antiviral(
     avs_itt_dat
)

n (\%)	Patients	Known	Missing	Shortness of breath day 28
Standard of care	73	67 (91.8)	6 (8.2)	35 (52.2)
Nafamostat	82	73 (89.0)	9 (11.0)	36 (49.3)
Overall	155	140 (90.3)	15 (9.7)	71 (50.7)

Table 2: Primary outcome by anti-viral intervention.

Age

Relationship age to shortness of breath

agedat <- fas_itt_dat %>%
  dplyr::count(out_sob, AgeAtEntry) %>% 
  spread(out_sob, n, fill = 0) %>% 
  mutate(
    n = `0` + `1` + `<NA>`,
    p = `1` / (`1` + `0`))
agemod <- glm(
  cbind(`1`, `0`) ~ AgeAtEntry, 
  data = agedat, 
  family = binomial())
agedat <- agedat %>%
  mutate(
    ypred = predict(agemod, newdata = agedat, type = "response")
  )
p1 <- ggplot(agedat, aes(AgeAtEntry, n)) +
  geom_col(colour = "grey40", fill = "grey40") +
  geom_vline(xintercept = 60, linetype = 2) +
  labs(y = "Number of\nparticipants",
       x = "Age at entry")
p2 <- ggplot(agedat, aes(AgeAtEntry, p)) +
    geom_point() +
    geom_vline(xintercept = 60, linetype = 2) +
    geom_line(aes(y = ypred)) +
    labs(y = "Proportion\nshortness of breath\nday 28", x = "Age at entry")
p <- p1 | p2
path <- file.path("outputs", "figures", "outcomes", "secondary")
ggsave(file.path(path, "7-6-age.pdf"), p, height = 2, width = 6)
p

Figure 1: Relationship (logistic regression linear in age) between age at entry and shortness of breath at day 28.

Sex

Relationship sex to outcome

tdat <- fas_itt_dat %>%
  dplyr::count(Sex, out_sob) %>%
  group_by(Sex) %>%
  spread(out_sob, n, fill = 0) %>%
  mutate(
    n = `1` + `0` + `<NA>`,
    p_1 = `1` / (`1` + `0`),
    p_miss = `<NA>` / (`1` + `0` + `<NA>`)
  )
p1 <- ggplot(tdat, aes(Sex, n)) +
  geom_col() +
    labs(
      y = "Number of\nparticipants", 
      x = "Sex")
p2 <- ggplot(tdat, aes(Sex, p_1)) +
  geom_point() +
    labs(
      y = "Proportion\nshortness of breath\nday 28", 
      x = "Sex")  +
  ylim(0, NA)
p <- p1 | p2
path <- file.path("outputs", "figures", "outcomes", "secondary")
ggsave(file.path(path, "7-6-sex.pdf"), p, height = 2, width = 6)
p

Figure 2: Shortness of breath by sex.

Oxygen

Relationship oxygen to outcome

tdat <- fas_itt_dat %>%
  dplyr::count(supp_oxy2 = factor(supp_oxy2, labels = c("Not required", "Required")), out_sob) %>%
  group_by(supp_oxy2) %>%
  spread(out_sob, n, fill = 0) %>%
  mutate(
    n = `1` + `0` + `<NA>`,
    p_1 = `1` / (`1` + `0`),
    p_miss = `<NA>` / (`1` + `0` + `<NA>`)
  )
p1 <- ggplot(tdat, aes(supp_oxy2, n)) +
  geom_col() +
    labs(
      y = "Number of\nparticipants", 
      x = "Supplemental oxygen")
p2 <- ggplot(tdat, aes(supp_oxy2, p_1)) +
  geom_point() +
    labs(
      y = "Proportion\nshortness of breath\nday 28", 
      x = "Supplemental oxygen")  +
  ylim(0, NA)
p <- p1 | p2
path <- file.path("outputs", "figures", "outcomes", "secondary")
ggsave(file.path(path, "7-6-oxygen.pdf"), p, height = 2, width = 6)
p

Figure 3: Shortness of breath by oxygen.

Country

Relationship country to outcome

tdat <- fas_itt_dat %>%
  dplyr::count(Country = factor(
    Country, 
    levels = c("IN", "AU", "NP", "NZ"),
    labels = c("India", "Australia", "Nepal", "New\nZealand")), out_sob) %>%
  group_by(Country) %>%
  spread(out_sob, n, fill = 0) %>%
  mutate(
    n = `1` + `0` + `<NA>`,
    p_1 = `1` / (`1` + `0`),
    p_miss = `<NA>` / (`1` + `0` + `<NA>`)
  )
p1 <- ggplot(tdat, aes(Country, n)) +
  geom_col() +
    labs(
      y = "Number of\nparticipants", 
      x = "Country of enrolment")
p2 <- ggplot(tdat, aes(Country, p_1)) +
  geom_point() +
    labs(
      y = "Proportion\nshortness of breath\nday 28", 
      x = "Country of enrolment")  +
  ylim(0, NA)
p <- p1 | p2
path <- file.path("outputs", "figures", "outcomes", "secondary")
ggsave(file.path(path, "7-6-country.pdf"), p, height = 2, width = 6)
p

Figure 4: Shortness of breath by country.

Site

Relationship site to outcome

tdat <- all_dat %>%
  filter_fas_itt() %>%
  dplyr::count(
    Country_lab = Country,
    Site_lab = fct_infreq(Location),
    Country = factor(PT_CountryName, levels = c("India", "Australia", "Nepal", "New Zealand")),
    Site = PT_LocationName,
    out_sob) %>%
  group_by(Country, Site) %>%
  spread(out_sob, n, fill = 0) %>%
  mutate(
    n = `1` + `0` + `<NA>`,
    p_1 = `1` / (`1` + `0`),
    p_miss = `<NA>` / (`1` + `0` + `<NA>`)
  ) %>%
  ungroup()
p1 <- ggplot(tdat, aes(Site_lab, n)) +
  facet_grid( ~ Country, scales = "free_x", space = "free_x") +
  geom_col() +
    labs(
      y = "Number of\nparticipants", 
      x = "") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        panel.border = element_rect(fill = NA))
p2 <- ggplot(tdat, aes(Site_lab, p_1)) +
  facet_grid( ~ Country, scales = "free_x", space = "free_x") +
  geom_point() +
    labs(
      y = "Proportion\nshortness of breath\nday 28", 
      x = "Site of enrolment") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        panel.border = element_rect(fill = NA))
p <- p1 / p2
path <- file.path("outputs", "figures", "outcomes", "secondary")
ggsave(file.path(path, "7-6-country-site.pdf"), p, height = 4, width = 6.25)
p

Figure 5: Shortness of breath at day 28 by site within country.

Calendar Time

Relationship calendar date to outcome

caldat <- all_dat %>% 
  filter_fas_itt() %>%
  dplyr::count(out_sob, yr = year(RandDate), mth = month(RandDate)) %>% 
  spread(out_sob, n, fill = 0) %>% 
  mutate(p = `1` / (`1` + `0`),
         n = `1` + `0` + `<NA>`)
p1 <- ggplot(caldat, aes(mth, p)) +
  facet_grid( ~ yr, drop = T, scales = "free_x", space = "free") +
    geom_point() +
    labs(
      y = "Proportion\nshortness of breath\nday 28", 
      x = "Calendar date (month of year)") +
  scale_x_continuous(breaks = 1:12) +
  ylim(0, NA)
p2 <- ggplot(caldat, aes(mth, n)) +
  facet_grid( ~ yr, drop = T, scales = "free_x", space = "free") +
    geom_col() +
    labs(
      y = "Number of\nparticipants", 
      x = "Calendar date (month of year)") +
  scale_x_continuous(breaks = 1:12)
p <- p2 | p1
path <- file.path("outputs", "figures", "outcomes", "secondary")
ggsave(file.path(path, "7-6-calendar-time.pdf"), p, height = 2, width = 6)
p

Figure 6: Relationship between calendar date and the primary outcome.

Analyses

FAS-ITT

Code

res <- fit_primary_model(fas_itt_nona_dat, logistic_site_epoch, ctr = contr.equalprior, outcome = "out_sob")
res$drws$OR <- res$drws$OR[-1]
names(res$drws$AOR) <- "Nafamostat"
names(res$drws$COR) <- intervention_labels2()$CAssignment[-(1:2)]
names(res$drws$OR) <- c("Ineligible aspirin", "Age \u2265 60", "Female", "Oxygen requirement", "Australia/New Zealand", "Nepal")

Odds ratio summary table

save_tex_table(
  odds_ratio_summary_table(c(res$drws$AOR, res$drws$COR, res$drws$OR), "latex"),
  "outcomes/secondary/7-6-primary-model-fas-itt-summary-table")
odds_ratio_summary_table(c(res$drws$AOR, res$drws$COR, res$drws$OR))

Parameter	Median	95% CrI	Mean (SD)	Pr(OR < 1)
Nafamostat	0.71	(0.30, 1.56)	0.76 (0.33)	0.81
Intermediate-dose	0.78	(0.55, 1.10)	0.79 (0.14)	0.92
Low-dose with aspirin	1.17	(0.75, 1.82)	1.20 (0.27)	0.24
Therapeutic-dose	0.97	(0.37, 2.39)	1.07 (0.53)	0.53
Ineligible aspirin	1.19	(0.40, 3.19)	1.34 (0.74)	0.37
Age ≥ 60	2.10	(1.50, 2.94)	2.14 (0.37)	0.00
Female	0.99	(0.72, 1.35)	1.00 (0.16)	0.54
Oxygen requirement	1.23	(0.86, 1.75)	1.25 (0.23)	0.14
Australia/New Zealand	2.67	(0.72, 9.63)	3.32 (2.48)	0.07
Nepal	0.48	(0.12, 2.36)	0.66 (0.74)	0.83

Posterior summaries for model parameters (fixed-effects), FAS-ITT.

Odds ratio summary for epoch and site

p <- plot_epoch_site_terms(
  res$drws$gamma_epoch,
  res$drws$gamma_site,
  factor(res$dat$region_by_site, 
         labels = c("India", "Australia\nNew Zealand", "Nepal"))
)
pth <- file.path("outputs", "figures", "outcomes", "secondary", "7-6-primary-model-epoch-site-terms-fas-itt.pdf")
ggsave(pth, p, width = 6, height = 4.5)
p

Summary of epoch and site posterior odds ratios.

Code

p <- plot_or_densities(c(res$drws$AOR, res$drws$COR)) +
  labs(x = "Odds ratio (log scale)", y = "Comparison")
pth <- file.path("outputs", "figures", "outcomes", "secondary", "7-6-primary-model-fas-itt-odds-ratio-densities.pdf")
ggsave(pth, p, width = 6, height = 2.5)
p

Posterior densities for treatment comparisons.

Diagnostics

Code

res$fit$diagnostic_summary()

$num_divergent
[1] 0 0 0 0 0 0 0 0

$num_max_treedepth
[1] 0 0 0 0 0 0 0 0

$ebfmi
[1] 0.7667 0.8095 0.8068 0.8222 0.7507 0.7393 0.8046 0.7920

Trace plots

Code

mcmc_trace(res$drws["beta"])

Code

mcmc_trace(res$drws["gamma_site"])

Code

mcmc_trace(res$drws["gamma_epoch"])

Posterior Predictive Checks

Code

y_ppc <- res$drws$y_ppc
ppc_dat <- bind_cols(fas_itt_nona_dat, tibble(y_ppc = y_ppc))

grp_ppc <- function(grp) {
  ppc_dat %>%
    group_by(grp = !!grp) %>%
    summarise(
      mean_y = mean(out_sob),
      rvar_mean_y_ppc = rvar_mean(y_ppc)
    )
}

plot_grp_ppc <- function(dat, lab = "") {
  dat %>%
    ggplot(aes(y = grp, xdist = rvar_mean_y_ppc)) +
    stat_interval(size = 2) +
    geom_point(aes(x = mean_y, y = 1:nrow(dat)), colour = "red", shape = 23) +
    labs(
      x = "Posterior predictive\nproportion", 
      y = lab)  
}

ppc_A <- grp_ppc(sym("AAssignment"))
ppc_C <- grp_ppc(sym("CAssignment"))
ppc_ctry <- grp_ppc(sym("Country"))
ppc_epoch <- grp_ppc(sym("epoch"))
ppc_site <- ppc_dat %>%
  group_by(Country = Country, grp = site) %>%
  summarise(
    mean_y = mean(out_sob),
    rvar_mean_y_ppc = rvar_mean(y_ppc)
  )

p0 <- plot_grp_ppc(ppc_A, "Anti-\nviral") + labs(x = "")
p1 <- plot_grp_ppc(ppc_C, "Anti-\ncoagulation") + labs(x = "")
p2 <- plot_grp_ppc(ppc_ctry, "Country") + labs(x = "")
p3 <- plot_grp_ppc(ppc_epoch, "Epoch") + labs(x = "")
p4 <- plot_grp_ppc(ppc_site %>% filter(Country == "IN"), "Sites\nIndia")  + labs(x = "")
p5 <- plot_grp_ppc(ppc_site %>% filter(Country == "AU"), "Sites\nAustralia") + labs(x = "")
p6 <- plot_grp_ppc(ppc_site %>% filter(Country == "NP"), "Sites\nNepal")
p7 <- plot_grp_ppc(ppc_site %>% filter(Country == "NZ"), "Sites\nNZ")
p <- ((p3 | (p0 | p1) / p2) / 
        ( (p4 | p5) / (p6 | p7) + 
            plot_layout(heights = c(3, 1)) ) ) +
  plot_layout(heights = c(1, 1.5), guides = "collect") &
  theme(legend.position = "bottom")
pth <- file.path("outputs", "figures", "outcomes", "secondary",
                 "7-6-primary-model-fas-itt-ppc.pdf")
ggsave(pth, p, width = 6, height = 5.5)
p

ACS-ITT

Fit primary model

res <- fit_primary_model(acs_itt_nona_dat, logistic_site_epoch, ctr = contr.equalprior, outcome = "out_sob")
res$drws$OR <- res$drws$OR[-1]
names(res$drws$AOR) <- "Nafamostat"
names(res$drws$COR) <- intervention_labels2()$CAssignment[-(1:2)]
names(res$drws$OR) <- c("Ineligible aspirin", "Age \u2265 60", "Female", "Oxygen requirement", "Australia/New Zealand", "Nepal")

Odds ratio summary table

save_tex_table(
  odds_ratio_summary_table(c(res$drws$AOR, res$drws$COR, res$drws$OR), "latex"),
  "outcomes/secondary/7-6-primary-model-acs-itt-summary-table")
odds_ratio_summary_table(c(res$drws$AOR, res$drws$COR, res$drws$OR))

Parameter	Median	95% CrI	Mean (SD)	Pr(OR < 1)
Nafamostat	0.85	(0.34, 2.07)	0.94 (0.45)	0.64
Intermediate-dose	0.77	(0.54, 1.09)	0.78 (0.14)	0.93
Low-dose with aspirin	1.17	(0.75, 1.81)	1.20 (0.27)	0.24
Therapeutic-dose	0.95	(0.36, 2.35)	1.06 (0.52)	0.54
Ineligible aspirin	1.23	(0.41, 3.34)	1.39 (0.77)	0.36
Age ≥ 60	2.04	(1.45, 2.90)	2.07 (0.37)	0.00
Female	1.02	(0.75, 1.40)	1.04 (0.16)	0.44
Oxygen requirement	1.28	(0.88, 1.85)	1.30 (0.25)	0.10
Australia/New Zealand	2.28	(0.61, 8.37)	2.84 (2.16)	0.11
Nepal	0.50	(0.12, 2.46)	0.69 (0.69)	0.82

Posterior summaries for model parameters (fixed-effects), FAS-ITT.

Odds ratio summary for epoch and site

p <- plot_epoch_site_terms(
  res$drws$gamma_epoch,
  res$drws$gamma_site,
  factor(res$dat$region_by_site, 
         labels = c("India", "Australia\nNew Zealand", "Nepal"))
)
pth <- file.path("outputs", "figures", "outcomes", "secondary", "7-6-primary-model-epoch-site-terms-acs-itt.pdf")
ggsave(pth, p, width = 6, height = 4.5)
p

Summary of epoch and site posterior odds ratios.

Code

p <- plot_or_densities(c(res$drws$AOR, res$drws$COR)) +
  labs(x = "Odds ratio (log scale)", y = "Comparison")
pth <- file.path("outputs", "figures", "outcomes", "secondary", "7-6-primary-model-acs-itt-odds-ratio-densities.pdf")
ggsave(pth, p, width = 6, height = 2.5)
p

Posterior densities for treatment comparisons.

Diagnostics

Code

res$fit$diagnostic_summary()

$num_divergent
[1] 0 0 0 0 0 0 0 0

$num_max_treedepth
[1] 0 0 0 0 0 0 0 0

$ebfmi
[1] 0.7511 0.7915 0.7607 0.7888 0.7786 0.7750 0.7824 0.8094

Trace plots

Code

mcmc_trace(res$drws["beta"])

Code

mcmc_trace(res$drws["gamma_site"])

Code

mcmc_trace(res$drws["gamma_epoch"])

AVS-ITT

Code

res <- fit_primary_model(
  avs_itt_nona_dat, 
  logistic_mod, 
  ctr = contr.equalprior,
  outcome = "out_sob", 
  vars = c("agegte60", "sexF", "supp_oxy2", "crp_tertile"), 
  beta_sd_var = c(2.5, 2.5, 2.5, 2.5, 2.5, 2.5)
)
res$drws$OR <- res$drws$OR[-1]
names(res$drws$AOR) <- "Nafamostat"
names(res$drws$COR) <- intervention_labels2()$CAssignment[-(1:2)]
names(res$drws$OR) <- c("Age \u2265 60", "Female", "Oxygen requirement", "CRP (2nd tertile)", "CRP (3rd tertile)", "CRP (unknown)")

Odds ratio summary table

save_tex_table(
  odds_ratio_summary_table(c(res$drws$AOR, res$drws$COR, res$drws$OR), "latex"),
  "outcomes/secondary/7-6-primary-model-avs-itt-summary-table")
odds_ratio_summary_table(c(res$drws$AOR, res$drws$COR, res$drws$OR))

Parameter	Median	95% CrI	Mean (SD)	Pr(OR < 1)
Nafamostat	0.76	(0.35, 1.62)	0.82 (0.33)	0.76
Intermediate-dose	0.78	(0.32, 1.88)	0.86 (0.41)	0.71
Low-dose with aspirin	0.58	(0.13, 2.38)	0.75 (0.61)	0.77
Therapeutic-dose	1.65	(0.51, 5.57)	2.00 (1.36)	0.20
Age ≥ 60	3.04	(1.33, 7.34)	3.37 (1.57)	0.00
Female	1.16	(0.55, 2.51)	1.26 (0.51)	0.35
Oxygen requirement	2.56	(1.09, 6.20)	2.83 (1.32)	0.01
CRP (2nd tertile)	0.83	(0.32, 2.18)	0.94 (0.49)	0.65
CRP (3rd tertile)	0.73	(0.28, 1.89)	0.82 (0.42)	0.74
CRP (unknown)	1.08	(0.30, 3.88)	1.34 (0.98)	0.45

Posterior summaries for model parameters (fixed-effects), AVS-ITT.

Code

p <- plot_or_densities(c(res$drws$AOR)) +
  labs(x = "Odds ratio (log scale)", y = "Comparison")
pth <- file.path("outputs", "figures", "outcomes", "secondary", "7-6-primary-model-avs-itt-odds-ratio-densities.pdf")
ggsave(pth, p, width = 6, height = 2.5)
p

Posterior densities for treatment comparisons.

Diagnostics

Code

res$fit$diagnostic_summary()

$num_divergent
[1] 0 0 0 0 0 0 0 0

$num_max_treedepth
[1] 0 0 0 0 0 0 0 0

$ebfmi
[1] 1.0006 0.9552 0.9447 0.9590 1.0144 0.9676 0.9890 0.9826

Trace plots

Code

mcmc_trace(res$drws["beta"])

End of script.